;;;; some random number generators
;;return a sampler, which is a thunk which returns a sample.

#!r6rs

(load "./monte-carlo.ss")

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; 0. Basic distributions
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; (a) uniform U(a,b)
(define (uniform a b)
  (lambda ()
    (+ a (random (- b a)))))

;; (b) normal N(mu,sigma)
(define (normal mu sigma)
  (define-values (e pi)
    (values (exp 1) (* (asin 1) 2)))
  (define (two-randoms)
    (let* ([u1 (random 1.0)]
           [u2 (random 1.0)]
           [aux1 (sqrt (* -2 (log u1 e)))]
           [aux2 (* 2 pi u2)])
      (map (lambda (proc)
             (+ mu (* sigma aux1 (proc aux2))))
           (list cos sin))))
  (let ([container (two-randoms)])
    (lambda ()
      (if (null? container) (set! container (two-randoms)))
      (let ([result (car container)])
          (set! container (cdr container))
          result)
      )))

;; (c) cosine Cosine(mu delta)
;; Cosine 提供了在紧致区域的单峰分布
;; pdf f(x) = pi/2/delta * cos(pi/delta* (x-mu))
;;inverse of cdf is F^-1(u) = mu + delta/pi * arcsin(2u-1)
(define (cosine mu delta);by sec-1 below
  (let* ([inv-pi (/ 1. 2 (asin 1.0))]
         [invF (lambda (u)
                 (+ mu (* delta inv-pi (asin (- (* 2 u) 1)))))])
    (lambda ()
      (invF (random 1.0)))))

;; (d) right/left-half cosine: rcosine(mu delta), lcosine(mu delta)
(define (rcosine mu delta)
  (let* ([inv-pi (/ 1. 2 (asin 1.0))]
         [invF (lambda (u) u)])
    (display "-------to be continued------")))



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; 1. Inverse transform method
;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;given a cdf F(X), then its inverse, G(U),U from U(0,1) can generate a rand-seq of X
(define (inverse-transform G)
  ;;G = F^-1
  (G (random 1.0)))

#|
;;example: integrate 2x^2 in (0,1) => result: 2/3
;;choose pdf f(x) = 2x in (0,1) => cdr F(x)=x^2 => inverse G(u) = sqrt(u)

(let* ([sampler (lambda ()
                  (inverse-transform
                   (lambda (u) (sqrt u))))]
       [int-strm (monte-carlo (lambda (x) x) sampler)]
       [ordinary-method (integral-with-errors (lambda (x) (* 2 x x)) 1)])
  (list (stream-ref int-strm 100000)
        (stream-ref ordinary-method 100000))
  )
|#


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; 2. Composition method
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; If you find the pdf f(X) = p1*f1+p2*f2+..., p1+p2+..=1,
;; and f1,f2,.. are also pdfs and easy to generate,
;; Then you can generate x in f1 in prob-p1, in f2 in pro-p2, ..

(define (composition . pss)
  ;;a ps is a pair:(p_i sampler_i)
  ;;change! pss to ((p1 . s1) (p1+p2 . s2) ..)
  (let ([aur 0.])
    (for-each (lambda (ps)
                (set! aur (+ aur (car ps)))
                (set-car! ps aur))
              pss))
  (lambda ()
    (let ([p (random 1.0)])
      (let loop ([pss pss])
        (if (< p (caar pss))
            ((cdar pss))
            (loop (cdr pss)))))))



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; 3 Acceptance-Rejection method
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;This method assumes that target pdf f(x) is controled by a proposal pdf g(x) by f(x)<=Cg(x)
;;the method is: 1. generate x by g; 2. generate y by unitform[0,Cg(x)]; 3. check y <= f(x) (yes,ok;no,reject)
;; note: C >=1, and this method hopes C is as small as possible
(define (acceptance-rejection f g gsampler C)
  (lambda ()
    (let loop ()
      (let* ([x (gsampler)]
             [y (* C (g x) (random 1.0))])
        (if (<= y (f x)) x
            (loop)))))
  )

;;example: integrate 2x^3 in [0,1]; h(x)=x; precise result: 2/3
;;target_pdf: f(x)=2x, in [0,1]; g(x)=1, C=2

#|

|#

(let ([newsampler (acceptance-rejection
                   (lambda (x) (* 2 x))
                   (lambda (x) 1)
                   (lambda () (random 1.0))
                   2)])
  (stream-ref (monte-carlo (lambda (x) x) newsampler)
              10000))


;; 4 Markov Chain Monte Carlo (MCMC)

(define (Matropolis-Hastings f Q Qsampler init);=> a sampler
  ;;f(x): target pdf(up to a constant); init: first item in the markov chain;
  ;;Q(x,y): proposal transition; Qsampler(x): sampler by Q(x, .)
  (let* ([rho (lambda (x y) (* (/ (f y) (f x))
                               (/ (Q y x) (Q x y))))]
         [alpha (lambda (x y) (min (rho x y) 1))])
    (let ([next-sample (lambda (current)
                         (let ([y ((Qsampler current))]
                               [u (random 1.0)])
                           (if (<= u (alpha current y))
                               y current)))]
          [current init])
      ;;过滤前面充分多项
      (do ((j 0 (+ j 1)))
          ((= j 10000))
        (set! current (next-sample current)))
      (lambda ()
        (set! current (next-sample current))
        (next-sample current)))))

;;test
(let* ([f (lambda (x) (sqrt (- 1 (expt x 2))))]
       [g (lambda (x) 1)]
       [Q (lambda (x y) (g y))]
       [Qsampler (lambda (x) (lambda () (random 1.0)))]
       [MH-sampler (Matropolis-Hastings f Q Qsampler (random 1.0))])
  (/ 4.0 (car (stream-ref (monte-carlo (lambda (x) (/ 1.0 (sqrt (- 1. (* x x)))))
                                       MH-sampler 'kappa)
                          1000000))))













